---
title: "Replication Analysis: From Criticism to Pragmatism? The Impact of NGEU Funding on Party-Based Euroscepticism in Italy"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 3
    toc-title: "Table of Contents"
    toc-location: left
    code-fold: true
editor: visual
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
rm(list = ls())

library(readxl)
library(tidyverse)
library(modelsummary)
library(marginaleffects)
library(here)
library(patchwork)
library(tinytable)
library(parameters)
library(forcats)
library(lme4)

theme_set(theme_bw())
options(scipen = 999)
options(marginaleffects_safe = FALSE)
```

```{r}
# Load the data
data <- read_csv(here("data_article.csv"))

data <- data |> 
  mutate(NRRP_Period_2020 = factor(
      NRRP_Period_2020,
      levels =c("Before NGEU", "After NGEU")
    ),
    NRRP_Period_2021 = factor(
      NRRP_Period_2021,
      levels = c("Before NGEU", "After NGEU")
    ),
    in_government = factor(
      in_government,
      levels = c("Not in Government", "In Government")
    ))

## Remove small parties and create data_small ####
data_small <- data %>%
  filter(!(
    party_cod %in% c(
      "Südtiroler Volkspartei",
      "More Europe",
      "Us Moderates",
      "Action", 
      "Italy Alive"
    )
  )) |>
  mutate(party_cod = factor(party_cod, levels = c("Left", "Democratic Party",
  "Five Star Movement", "Go Italy", "Brothers of Italy", "League")),
  party_cod = factor(party_cod, labels = c("SI", "PD", "M5S", "FI", "FdI", "Lega")))
```

# Descriptive statistics

## Table 1. Type of parliamentary statements, by party

```{r}
## all parties no policy/polity####
dat_summary <- data |> 
  ungroup() |>
  group_by(party_cod) |>
  summarise(n = n(),
            support = sum(support_cat, na.rm = TRUE),
            critique = sum(critique_cat, na.rm = TRUE),
            alternative = sum(alternative_cat, na.rm = TRUE))

long_data <- dat_summary %>%
  select(party_cod, n, support:alternative) |>
  pivot_longer(support:alternative,
               names_to = "facet_type",
               values_to = "count") |>
  mutate(facet_type = factor(
    facet_type,
    levels = c("support", "critique", "alternative"),
    labels = c("Support", "Critique", "Alternative"),
    ordered = TRUE
  )) 

long_data |> 
  pivot_wider(
    names_from = c(facet_type),
    values_from = count,
    values_fill = 0
  ) |> 
  mutate(
    Support_pct = round((Support / n) * 100, 1),
    Critique_pct = round((Critique / n) * 100, 1),
    Alternative_pct = round((Alternative / n) * 100, 1)
  ) |> 
  arrange(desc(n)) %>%
  mutate(
    Support_combined = paste0(Support, "(", round((Support / n) * 100, 1), "%)"),
    Critique_combined = paste0(Critique, "(", round((Critique / n) * 100, 1), "%)"),
    Alternative_combined = paste0(Alternative, "(", round((Alternative / n) * 100, 1), "%)")
  ) %>%
  select(party_cod, n, Support_combined, Critique_combined, Alternative_combined) |> 
  tt() |>
  setNames(c("Party", "Statements", "Support", "Critique", "Alternative"))
```

# Model estimation

```{r}
explanatory_formula_noint <- "party_cod + NRRP_Period_2020 + in_government"
explanatory_formula_int <- "party_cod * NRRP_Period_2020 + in_government"
explanatory_formula_govparty <- "NRRP_Period_2020 * in_government + party_cod"
explanatory_formula_sup <- "party_cod + eusupport_index + in_government"
explanatory_formula_sup_int <- "party_cod * eusupport_index + in_government"
explanatory_formula_2021 <- "party_cod * NRRP_Period_2021 + in_government"
explanatory_formula_noint_2021 <- "party_cod + NRRP_Period_2021 + in_government"

model_run <- function(explanatory_form) {
  model_support <- glm(as.formula(paste("support_cat ~", explanatory_form)), family = binomial, data = data_small)
  model_critique <- glm(as.formula(paste("critique_cat ~", explanatory_form)), family = binomial, data = data_small)
  model_alternative <- glm(as.formula(paste("alternative_cat ~", explanatory_form)), family = binomial, data = data_small)
  return(list(support = model_support, critique = model_critique, alternative = model_alternative))
}

model_noint <- model_run(explanatory_formula_noint)
model_int <- model_run(explanatory_formula_int)
model_govparty <- model_run(explanatory_formula_govparty)
model_sup <- model_run(explanatory_formula_sup)
model_sup_int <- model_run(explanatory_formula_sup_int)
model_2021 <- model_run(explanatory_formula_2021)
model_noint_2021 <- model_run(explanatory_formula_noint_2021)
```

## Table 2

```{r}

models <-   list(
  "Support" = model_noint[[1]],
  "Criticism" = model_noint[[2]],
  "Alternative" = model_noint[[3]],
  "Support (with interaction)" = model_int[[1]],
  "Criticism (with interaction)" = model_int[[2]],
  "Alternative (with interaction)" = model_int[[3]]
)

cm <- c('NRRP_Period_2020After NGEU' = 'NGEU',
        'in_governmentIn Government' = 'In government',
        'party_codPD'    = 'PD',
        'party_codM5S' = 'M5S',
        'party_codFI'  = 'FI',
        'party_codFdI'  = 'FdI',
        'party_codLega' = 'League',
        'party_codPD:NRRP_Period_2020After NGEU' = 'PD × NGEU',
        'party_codM5S:NRRP_Period_2020After NGEU' = 'M5S × NGEU',
        'party_codFI:NRRP_Period_2020After NGEU' = 'FI × NGEU',
        'party_codFdI:NRRP_Period_2020After NGEU' = 'FdI × NGEU',
        'party_codLega:NRRP_Period_2020After NGEU' = 'League × NGEU')

modelsummary(
  models,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  vcov = ~date_chr,
  coef_omit = "Intercept",
  coef_map = cm,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```

### H1. Support

```{r}
# 1. Calculate the average predicted probability of making a support statement before and after NGEU (using clustered SEs by date).
avg_predictions(model_noint[[1]], by = "NRRP_Period_2020", vcov = ~date_chr)
avg_predictions(model_noint[[1]], by = "NRRP_Period_2020", vcov = ~date_chr, hypothesis = "b2-b1 = 0")
pred_support <- plot_predictions(model_noint[[1]], by = "NRRP_Period_2020", vcov = ~date_chr) + labs(title="Support statements", x = "", y = "")
```

```{r}
avg_slopes(model_noint[[1]], variables = "NRRP_Period_2020", vcov = ~date_chr)
```

### H1. Critique

```{r}
avg_predictions(model_noint[[2]], by = "NRRP_Period_2020", vcov = ~date_chr)
avg_predictions(model_noint[[2]], by = "NRRP_Period_2020", vcov = ~date_chr, hypothesis = "b2-b1 = 0")
pred_critique <- plot_predictions(model_noint[[2]], by = "NRRP_Period_2020", vcov = ~date_chr) + labs(title="Critique statements", x = "", y = "")
```

```{r}
avg_slopes(model_noint[[2]], variables = "NRRP_Period_2020", vcov = ~date_chr)
```

### H1. Alternative

```{r}
avg_predictions(model_noint[[3]], by = "NRRP_Period_2020", vcov = ~date_chr)
avg_predictions(model_noint[[3]], by = "NRRP_Period_2020", vcov = ~date_chr, hypothesis = "b2-b1 = 0")
pred_alternative <- plot_predictions(model_noint[[3]], by = "NRRP_Period_2020", vcov = ~date_chr) + labs(title="Alternative statements", x = "", y = "")
```

```{r}
avg_slopes(model_noint[[3]], variables = "NRRP_Period_2020", vcov = ~date_chr)
```

## Figure 1. Average predicted probability of making a statement about the EU before and after NGEU

```{r}
pred_support + pred_critique + pred_alternative
```

### H2. Support

```{r}
avg_predictions(model_int[[1]], by = c("party_cod", "NRRP_Period_2020"), vcov = ~date_chr)
```

```{r}
avg_slopes(model_int[[1]], variables = "NRRP_Period_2020", by = "party_cod", vcov = ~date_chr)
```

### H2. Critique

```{r}
avg_predictions(model_int[[2]], by = c("party_cod", "NRRP_Period_2020"), vcov = ~date_chr)
```

```{r}
avg_slopes(model_int[[2]], variables = "NRRP_Period_2020", by = "party_cod", vcov = ~date_chr)
```

### H2. Alternative

```{r}
avg_predictions(model_int[[3]], by = c("party_cod", "NRRP_Period_2020"), vcov = ~date_chr)
```

```{r}
avg_slopes(model_int[[3]], variables = "NRRP_Period_2020", by = "party_cod", vcov = ~date_chr)
```

## Figures 2-4 Average predicted probability of making a support/critical/alternative statement about the EU before and after NGEU by party

```{r}
analyze_margin <- function(model) {
  library(ggplot2)
  library(marginaleffects)
  
  plot_predictions(
    model,
    by = c("party_cod", "NRRP_Period_2020"),
    vcov = ~date_chr
  ) +
    labs(
      x = "Party",
      y = "Average predicted probability",
      colour = "", # Optional if `by` is used for shapes
    ) +
    theme_bw() +
    scale_color_manual(values = c("black", "gray50", "gray80")) + # Black and white colors
    scale_fill_manual(values = c("black", "gray50", "gray80")) + # Ensure consistency
    theme(
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.text = element_text(size = 8),
      axis.title = element_text(size = 12),
      axis.text = element_text(size = 10)
    )
}

analyze_margin(model = model_int[[1]])
ggsave("model_support_margin.png", dpi = 300)

analyze_margin(model = model_int[[2]])
ggsave("model_critique_margin.png", dpi = 300)

analyze_margin(model = model_int[[3]])
ggsave("model_alternative_margin.png", dpi = 300)
```

# ANNEXES

## Table A3-A8. Combinations of statement types by party

```{r}
# Create all 2x2x2 combinations for display
# All possible combinations (0/1 for each of the three variables)
combinations <- expand.grid(
  support_cat = 0:1,
  critique_cat = 0:1,
  alternative_cat = 0:1
)

parties <- unique(as.character(data_small$party_cod))

party_tables <- map(
  parties,
  function(parties) {
    party_counts <- data_small %>%
      filter(party_cod == parties) %>%
      group_by(support_cat, critique_cat, alternative_cat) %>%
      summarise(n = n(), .groups = "drop") %>%
      right_join(combinations, by = c("support_cat", "critique_cat", "alternative_cat")) %>%
      mutate(n = replace_na(n, 0)) %>%
      mutate(
        Combination = paste0(
          "Support=", support_cat,
          " + Critique=", critique_cat,
          " + Alternative=", alternative_cat
        )
      ) %>%
      select(Combination, n) |> 
   knitr::kable(caption = paste0("Combinations of statement types for ", parties))
  }
)
names(party_tables) <- parties

walk(party_tables, print)
```

## Table A9. Model testing the impact of public support for the EU (data drawn from eurobarometer)

```{r}
models_sup <-   list(
  "Support" = model_sup[[1]],
  "Criticism" = model_sup[[2]],
  "Alternative" = model_sup[[3]],
  "Support (with interaction)" = model_sup_int[[1]],
  "Criticism (with interaction)" = model_sup_int[[2]],
  "Alternative (with interaction)" = model_sup_int[[3]]
)

cm_sup <- c(        'in_governmentIn Government' = 'In government',
        'eusupport_index' = 'EU support',
        'party_codPD'    = 'PD',
        'party_codM5S' = 'M5S',
        'party_codFI'  = 'FI',
        'party_codFdI'  = 'FdI',
        'party_codLega' = 'Lega',
        'party_codPD:eusupport_index' = 'PD × EU support',
        'party_codM5S:eusupport_index' = 'M5S × EU support',
        'party_codFI:eusupport_index' = 'FI × EU support',
        'party_codFdI:eusupport_index' = 'FdI × EU support',
        'party_codLega:eusupport_index' = 'Lega × EU support')

modelsummary(
  models_sup,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  vcov = ~date_chr,
  coef_omit = "Intercept",
  coef_map = cm_sup,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```

```{r}
avg_slopes(model_sup[[1]], variables = "eusupport_index", vcov = ~date_chr)
avg_slopes(model_sup[[2]], variables = "eusupport_index", vcov = ~date_chr)
avg_slopes(model_sup[[3]], variables = "eusupport_index", vcov = ~date_chr)

```

## Table A10. Focus on FdI in opposition.

```{r}
data_small %>%
  filter(government != "Meloni", party_cod == "FdI") %>% 
  group_by(NRRP_Period_2020) %>%
  summarise(
    Total = n(),
    Support = paste0(sum(support_cat, na.rm = TRUE), 
                     " (", round(100 * sum(support_cat, na.rm = TRUE) / Total, 1), "%)"),
    Critique = paste0(sum(critique_cat, na.rm = TRUE), 
                      " (", round(100 * sum(critique_cat, na.rm = TRUE) / Total, 1), "%)"),
    Alternative = paste0(sum(alternative_cat, na.rm = TRUE), 
                         " (", round(100 * sum(alternative_cat, na.rm = TRUE) / Total, 1), "%)"),
    .groups = "drop"
  ) %>%
  arrange(NRRP_Period_2020) %>%
  select(NRRP_Period_2020, Support, Critique, Alternative, Total) %>%  # <- Total at the end
  tt(notes = "Removed the statements during the Meloni government")
```

## Table A11. Focus on Lega in government

```{r}

data_small %>%
  filter(government != "Conte II" & party_cod == "Lega") %>% 
  group_by(NRRP_Period_2020) %>%
  summarise(
    Total = n(),
    Support = paste0(sum(support_cat, na.rm = TRUE), 
                     " (", round(100 * sum(support_cat, na.rm = TRUE) / Total, 1), "%)"),
    Critique = paste0(sum(critique_cat, na.rm = TRUE), 
                      " (", round(100 * sum(critique_cat, na.rm = TRUE) / Total, 1), "%)"),
    Alternative = paste0(sum(alternative_cat, na.rm = TRUE), 
                         " (", round(100 * sum(alternative_cat, na.rm = TRUE) / Total, 1), "%)"),
    .groups = "drop"
  ) %>%
  arrange(NRRP_Period_2020) %>%
  select(NRRP_Period_2020, Support, Critique, Alternative, Total) %>%
  tt(notes = "Removed the statements during the Conte II government")
```

## Table A12. Focus on M5S in government

```{r}
data_small %>%
  filter(government != "Meloni" & party_cod == "M5S") %>% 
  group_by(NRRP_Period_2020) %>%
  summarise(
    Total = n(),
    Support = paste0(sum(support_cat, na.rm = TRUE), 
                     " (", round(100 * sum(support_cat, na.rm = TRUE) / Total, 1), "%)"),
    Critique = paste0(sum(critique_cat, na.rm = TRUE), 
                      " (", round(100 * sum(critique_cat, na.rm = TRUE) / Total, 1), "%)"),
    Alternative = paste0(sum(alternative_cat, na.rm = TRUE), 
                         " (", round(100 * sum(alternative_cat, na.rm = TRUE) / Total, 1), "%)"),
    .groups = "drop"
  ) %>%
  arrange(NRRP_Period_2020) %>%
  select(NRRP_Period_2020, Support, Critique, Alternative, Total) %>% 
 tt(notes = "Removed the statements during the Meloni government")
```

## Table A13. Focus on CONTE II

```{r}

data_small %>%
  filter(government == "Conte II") %>% # Filter for Conte II government
  group_by(party_cod, NRRP_Period_2020) %>%
  summarise(
    Total = n(),
    Support = paste0(sum(support_cat, na.rm = TRUE), 
                     " (", round(100 * sum(support_cat, na.rm = TRUE) / Total, 1), "%)"),
    Critique = paste0(sum(critique_cat, na.rm = TRUE), 
                      " (", round(100 * sum(critique_cat, na.rm = TRUE) / Total, 1), "%)"),
    Alternative = paste0(sum(alternative_cat, na.rm = TRUE), 
                         " (", round(100 * sum(alternative_cat, na.rm = TRUE) / Total, 1), "%)"),
    .groups = "drop"
  ) %>%
  arrange(party_cod, NRRP_Period_2020)|> 
  select(party_cod, NRRP_Period_2020, Support, Critique, Alternative, Total) %>%  # <- Total at the end
  tt()
```

## Table A14. Model testing the interaction of government status and NGEU - Binomial logistic model with fixed effects for party

```{r}
models_govparty <-   list(
  "Support" = model_govparty[[1]],
  "Criticism" = model_govparty[[2]],
  "Alternative" = model_govparty[[3]]
)

cm_govparty <- c(        'NRRP_Period_2020After NGEU' = 'NGEU',
        'in_governmentIn Government' = 'In government',
        'party_codPD'    = 'PD',
        'party_codM5S' = 'M5S',
        'party_codFI'  = 'FI',
        'party_codFdI'  = 'FdI',
        'party_codLega' = 'Lega',

        'NRRP_Period_2020After NGEU:in_governmentIn Government'
        = 'NGEU × In government')

modelsummary(
  models_govparty,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  vcov = ~date_chr,
  coef_omit = "Intercept",
  coef_map = cm_govparty,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```

```{r}
library(marginaleffects)

analyze_govparty <- function(model) {
  # Calculate average comparisons
  comp1 <- avg_comparisons(
    model,
    variables = "in_government",
    by = "NRRP_Period_2020",
    vcov = ~date_chr
  )
  
  # Perform hypothesis testing
  comp2 <- marginaleffects::hypotheses(
    comp1,
    hypothesis = "b2 - b1 = 0"
  )
  
  # Return both results
  return(list(comp1 = comp1, comp2 = comp2))
}

analyze_govparty(model = model_govparty[[1]])
analyze_govparty(model = model_govparty[[2]])
analyze_govparty(model = model_govparty[[3]])
```

## Table A15. Model testing the interaction of government status and NGEU - Binomial logistic mixed-effects models with random intercepts for party–cabinet

```{r}
# ---- Prep: factors + party–cabinet ----
data_model <- data_small %>%
  mutate(
    in_government = fct_relevel(as.factor(in_government),
                                "Not in Government", "In Government"),
    NRRP_Period_2020 = fct_relevel(as.factor(NRRP_Period_2020),
                                   "Before NGEU", "After NGEU"),
    party_cabinet = fct_drop(interaction(party_cod, government, drop = TRUE))
  )

# Outcomes to run (names = column names; values = pretty labels)
outcomes <- c(
  support_cat    = "Support",
  critique_cat   = "Critique",
  alternative_cat= "Alternative"
)

# ---- Helper: fit a GLMM for a given outcome ----
fit_glmm <- function(outcome) {
  glmer(
    as.formula(paste0(outcome, " ~ in_government * NRRP_Period_2020 + (1 | party_cabinet)")),
    data = data_model,
    family = binomial(link = "logit"),
    control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
  )
}

# ---- Helper: predictions (population level) with robust column names ----
pred_2x2 <- function(model, outcome_label) {
  predictions(
    model,
    newdata = datagrid(
      in_government = c("Not in Government", "In Government"),
      NRRP_Period_2020 = c("Before NGEU", "After NGEU")
    ),
    re.form = NA,
    type = "response"
  ) %>%
    # harmonize possible column name variants from marginaleffects
    rename(
      estimate  = any_of(c("Estimate", "estimate", "predicted")),
      conf.low  = any_of(c("2.5 %", "conf.low")),
      conf.high = any_of(c("97.5 %", "conf.high"))
    ) %>%
    mutate(
      in_government = factor(in_government, levels = c("Not in Government", "In Government")),
      NRRP_Period_2020 = factor(NRRP_Period_2020, levels = c("Before NGEU", "After NGEU")),
      Outcome = outcome_label
    )
}

# ---- Fit all models + collect predictions ----
models_govparty2 <- map(names(outcomes), fit_glmm) |> set_names(outcomes)  # names become pretty labels

cm_govparty2 <- c(        'NRRP_Period_2020After NGEU' = 'NGEU',
        'in_governmentIn Government' = 'In government',
        'in_governmentIn Government:NRRP_Period_2020After NGEU'
        = 'NGEU × In government')

modelsummary(
  models_govparty2,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  coef_omit = "Intercept",
  coef_map = cm_govparty2,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```

## Table A16. Models Excluding May–December 2020

```{r}

data_no2020 <- data_small |> 
  filter(date_chr < as.Date("2020-05-01") | date_chr > as.Date("2020-12-31"))

model_run_no2020 <- function(explanatory_form) {
  model_support <- glm(as.formula(paste("support_cat ~", explanatory_form)), family = binomial, data = data_no2020)
  model_critique <- glm(as.formula(paste("critique_cat ~", explanatory_form)), family = binomial, data = data_no2020)
  model_alternative <- glm(as.formula(paste("alternative_cat ~", explanatory_form)), family = binomial, data = data_no2020)
  return(list(support = model_support, critique = model_critique, alternative = model_alternative))
}

model_noint_no2020 <- model_run_no2020(explanatory_formula_noint)
model_int_no2020 <- model_run_no2020(explanatory_formula_int)
model_govparty_no2020 <- model_run_no2020(explanatory_formula_govparty)
model_sup_no2020 <- model_run_no2020(explanatory_formula_sup)
model_sup_int_no2020 <- model_run_no2020(explanatory_formula_sup_int)
model_2021_no2020 <- model_run_no2020(explanatory_formula_2021)
model_noint_2021_no2020 <- model_run_no2020(explanatory_formula_noint_2021)

models_no2020 <-   list(
  "Support" = model_noint_no2020[[1]],
  "Criticism" = model_noint_no2020[[2]],
  "Alternative" = model_noint_no2020[[3]],
  "Support (with interaction)" = model_int_no2020[[1]],
  "Criticism (with interaction)" = model_int_no2020[[2]],
  "Alternative (with interaction)" = model_int_no2020[[3]]
)

cm_no2020 <- c('NRRP_Period_2020After NGEU' = 'NGEU',
        'in_governmentIn Government' = 'In government',
        'party_codPD'    = 'PD',
        'party_codM5S' = 'M5S',
        'party_codFI'  = 'FI',
        'party_codFdI'  = 'FdI',
        'party_codLega' = 'League',
        'party_codPD:NRRP_Period_2020After NGEU' = 'PD × NGEU',
        'party_codM5S:NRRP_Period_2020After NGEU' = 'M5S × NGEU',
        'party_codFI:NRRP_Period_2020After NGEU' = 'FI × NGEU',
        'party_codFdI:NRRP_Period_2020After NGEU' = 'FdI × NGEU',
        'party_codLega:NRRP_Period_2020After NGEU' = 'League × NGEU')

modelsummary(
  models_no2020,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  vcov = ~date_chr,
  coef_omit = "Intercept",
  coef_map = cm,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```

## Table A17. Model with alternative temporal cutoff point (2021-08-13 European Commission disburses €24.9 billion in pre-financing to Italy)

```{r}
models_2021 <-   list(
  "Support" = model_noint_2021[[1]],
  "Criticism" = model_noint_2021[[2]],
  "Alternative" = model_noint_2021[[3]],
  "Support (with interaction)" = model_2021[[1]],
  "Criticism (with interaction)" = model_2021[[2]],
  "Alternative (with interaction)" = model_2021[[3]]
)

cm2021 <- c('NRRP_Period_2021After NGEU' = 'NGEU',
        'in_governmentIn Government' = 'In government',
        'party_codPD'    = 'PD',
        'party_codM5S' = 'M5S',
        'party_codFI'  = 'FI',
        'party_codFdI'  = 'FdI',
        'party_codLega' = 'Lega',
        'party_codPD:NRRP_Period_2021After NGEU' = 'PD × NGEU',
        'party_codM5S:NRRP_Period_2021After NGEU' = 'M5S × NGEU',
        'party_codFI:NRRP_Period_2021After NGEU' = 'FI × NGEU',
        'party_codFdI:NRRP_Period_2021After NGEU' = 'FdI × NGEU',
        'party_codLega:NRRP_Period_2021After NGEU' = 'Lega × NGEU')

modelsummary(
  models_2021,
  estimate = "{estimate} ({std.error}){stars}",
  fmt = 2,
  statistic = NULL,
  vcov = ~date_chr,
  coef_omit = "Intercept",
  coef_map = cm2021,
  title = NULL,
  gof_omit = c("AIC|BIC|Std.Errors|RMSE"))
```
